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In  this  report,  we  describe  a  fast  deconvolution  approach  for  color  images  that  combines  a  sparse  regularization  cost 
on  the  magnitudes  of  gradients  with  constraints  on  their  direction  in  color  space.  We  form  these  color  constraints  in  a 
way  that  allows  retaining  the  computationally-efficient  optimization  strategy  introduced  in  recent  deconvolution 
methods  based  on  \emph  {half-quadratic  splitting}.  The  proposed  algorithm  is  capable  of  handling  a  different  blur 
kernel  in  each  color  channel,  and  is  used  for  per-layer  deconvolution  in  our  paper:  Depth  and  Deblurring  from  a 
Spectrally- varying  Depth-of-Field.  A  MATLAB  implementation  of  this  method  takes  roughly  20  seconds  to 
deconvolve  a  three-channel  1544x1028  color  image,  on  a  Linux-based  Intel  1-3  2.1GHz  machine. 
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Abstract — In  this  report,  we  describe  a  fast  deconvolution  approach 
for  color  images  that  combines  a  sparse  regularization  cost  on  the 
magnitudes  of  gradients  with  constraints  on  their  direction  in  color 
space.  We  form  these  color  constraints  in  a  way  that  allows  retaining 
the  computationally-efficient  optimization  strategy  introduced  in  recent 
deconvolution  methods  based  on  half-quadratic  splitting.  The  proposed 
algorithm  is  capable  of  handling  a  different  blur  kernel  in  each  color 
channel,  and  is  used  for  per-layer  deconvolution  in  our  paper:  “Depth 
and  Deblurring  from  a  Spectrally-varying  Depth-of-Field”  [1].  A  MATLAB 
implementation  of  this  method  is  available  at  http://vision.seas.harvard. 
edu/ccap,  and  takes  roughly  20  seconds  to  deconvolve  a  three-channel 
1544  x  1028  color  image,  on  a  Linux-based  Intel  1-3  2.1GHz  machine. 

1  Introduction 

Non-blind  deconvolution  refers  to  the  recovery  of  a  sharp 
image  x(n)  from  a  blurred  and  noisy  observation  y(n): 

y(n)  =  (x  *  k)(n)  +  e(n),  (1) 

where  the  blur  kernel  k  is  known,  and  e(n)  is  typically  assumed 
to  be  white  Gaussian  noise.  A  regularized  estimation  approach 
recovers  the  sharp  image  x(n)  as 

x(n)  =  argmin  ^  V  [y(n)  —  {x  *  k)(n)]2  +  $(x),  (2) 

x  2 

n 

where  the  first  term  measures  fidelity  to  the  observation,  <t>(x) 
is  a  regularization  cost  based  on  statistical  properties  of  sharp 
natural  images,  and  p  is  the  relative  weight  between  the  two. 

A  common  choice  for  $(#)  is  a  smoothness  cost  that  dis¬ 
courages  discontinuities  in  the  sharp  image  x  by  penalizing 
gradient  magnitudes  [2],  [3],  [4]  as 

$(i)=  (3) 

VeG  n 

where  a  >  0,  Q  is  a  set  of  gradient  filters,  and  x^(n)  =  (V  * 
x){n).  For  a  =  2,  the  solution  for  (2)  has  a  closed  form  that 
can  be  computed  efficiently  in  the  Fourier  domain.  However, 
such  a  squared  cost  prefers  to  distribute  gradients  equally  over 
the  recovered  image  [2],  and  based  on  the  choice  of  /i,  yields 
images  that  are  either  over-smoothed,  or  have  ringing  artifacts. 

It  is  generally  understood  that  a  choice  of  a  E  (0, 1]  corre¬ 
sponds  better  to  the  distribution  of  gradients  in  natural  images, 
and  yields  solutions  with  sharper  edges  where  gradients  are 
concentrated  at  a  sparse  set  of  pixels.  Unfortunately,  these 
values  of  a  do  not  admit  closed-form  solutions  and  require  an 
iterative  approach  for  deconvolution.  Levin  et  al.  [2]  propose 
such  an  approach  based  on  iterative  re-weighted  least  squares 
(IRLS),  where  at  each  iteration  the  regularizer  is  approximated 
with  a  weighted  squared  cost  with  weights  based  on  the 
current  estimate  of  x.  Unfortunately  with  spatially-varying 
weights,  deconvolution  with  the  approximated  squared  cost 


can  not  be  carried  out  in  the  Fourier  domain,  and  the  solution 
at  each  iteration  must  in  turn  be  computed  iteratively  using 
the  conjugate-gradient  method. 

Recent  work  [3],  [4]  demonstrates  that  one  can  solve  (2)  for 
a  <  2  with  greater  efficiency  using  an  optimization  approach 
based  on  half-quadratic  splitting  (which  we  describe  in  de¬ 
tail  in  Sec.  2.1).  These  algorithms  involve  iterating  between 
per-pixel  shrinkage  operations  on  gradients,  and  deconvolu¬ 
tion  in  the  Fourier  domain.  A  comparison  in  [4]  shows  that 
this  approach  is  about  400  times  faster  than  IRLS  on  large 
"megapixel"  images  (i.e.,  those  with  a  resolution  greater  than 
1024  x  1024).  In  this  work,  we  extend  the  half-quadratic  split¬ 
ting  approach  to  include  color  constraints  on  gradients  in  the 
sharp  image  x  in  a  way  that  yields  improved  results,  with 
roughly  the  same  computational  cost. 

1.1  Color-based  Regularization 

The  regularization  schemes  discussed  above  are  defined  exclu¬ 
sively  in  terms  of  single-channel  image  statistics,  and  ignore 
color  information.  Recently,  Joshi  et  al.  [5]  demonstrated  that 
higher  quality  results  can  be  obtained  by  using  color  statistics 
during  deconvolution,  even  for  the  case  when  the  blur  kernel 
is  spectrally  uniform.  Color  models  are  even  more  beneficial 
in  the  context  of  deconvolution  in  our  work  in  [1],  where 
each  color  channel  is  blurred  with  a  different  kernel.  In  this 
case,  observed  artifacts  from  independent  deconvolution  of 
color  channels  correspond  to  different  spatial  frequencies  in 
the  different  channels,  and  can  lead  to  the  creation  of  spurious 
chromaticities  in  the  estimated  image.  Conversely,  by  incor¬ 
porating  color  statistics,  one  can  exploit  the  fact  the  observed 
image  contain  complementary  spatial  frequency  information  in 
the  different  channels. 

The  color-based  regularization  introduced  in  [5]  is  defined  in 
terms  of  pixel-domain  statistics  in  local  neighborhoods,  and  the 
corresponding  estimation  algorithm  is  developed  in  the  IRLS 
framework.  In  this  work,  we  aim  to  incorporate  color  statistics 
in  a  way  that  admits  use  of  the  fast  deconvolution  framework 
of  [3],  [4].  We  reason  that  with  a  high  weight  on  the  regularizer 
during  estimation,  we  recover  an  image  Xs  that  is  excessively 
smoothed  but  has  fewer  artifacts  due  to  noise,  and  can  be 
considered  to  be  a  blurred  version  of  the  true  sharp  image 
A.  We  use  a  squared-cost  with  a  low  value  of  /i  to  estimate 
this  over-regularized  estimate  Xs  in  closed-form. 

Then,  we  draw  on  intuition  from  the  spatio-spectral  model 
in  [1]  that  gradients  in  local  neighborhoods  have  the  same 
color  direction,  i.e.,  they  lie  on  the  same  line  in  color  space. 
This  means  that  gradients  in  Xs/  being  local  averages  of  the 
true  gradients,  are  likely  to  lie  on  the  same  color  line  as 
those  in  X,  albeit  with  attenuated  magnitudes.  Our  strategy 
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is  to  use  the  gradient  colors  from  Xs  as  additional  constraints 
while  estimating  X  with  a  traditional  sparse  regularization  cost 
on  gradient  magnitudes.  We  find  that  this  simple  approach 
achieves  performance  comparable  to  that  of  the  color  IRLS- 
based  method  in  [5],  at  a  fraction  of  the  computational  cost. 

2  Proposed  Method 

We  consider  the  general  case  where  each  color  channel  y^(n) 
of  Y(n),i  G  {R,G,B}  is  blurred  with  a  different  kernel  . 
We  first  generate  an  over-regularized  estimate  Xs  of  X  with  a 
squared  cost  and  /xs  =  0.5  as: 

=  argmin  ^  |y^(n)  —  (x  *  fc^)(n)j 

n 

+  Y  Y  ia;v(n)i2  ’  (4) 

VeQ  n 

where  the  set  Q  contains  two  gradient  filters  corresponding 
to  vertical  and  horizontal  finite-differences.  This  can  be  solved 
directly  in  the  Fourier  domain  as: 

1 

Evl^[v]|2  +  ^[fc«]l2J  ’  u 

where  T  and  T~x  refer  to  the  2D  Fourier  transform  and  its 
inverse.  Note  that  this  Fourier-domain  formulation  is  based 
on  a  periodic-extension  assumption  on  the  image.  To  prevent 
ringing  from  sharp  discontinuities  between  pixels  at  opposite 
boundaries,  we  extend  the  observed  image  Y (n)  by  padding 
it  with  values  that  smoothly  blend  these  discontinuities.  These 
extended  regions  are  discarded  after  deconvolution  to  yield  a 
result  with  the  original  resolution. 

Xs  (n)  serves  as  the  reference  for  gradient  colors  in  the  final 
estimate  X.  For  notational  convenience,  we  define  Xv(n)  as  a 
color-gradient  vector  of  X  (n)  as: 

Xv(n)  =  [z^}(n),a4G}(n),z4S}(n)]  ,  (6) 

and  the  gradient  color  unit-vectors  Vv(n)  from  Xs  as 

Vv(n)  =  XaV(n)/||XaV(n)||.  (7) 

The  final  sharp  image  is  then  recovered  as 

X  =  argmin  ^  jy^(n)  —  {x^  *  fc^)(n)j 

n,i 

+  yyi(Jv(n),yv(n))|“,  (8) 

VeQ  n 

with  the  constraint  that  X^(n)  oc  Vv(n),  Vn.  We  use  a  =  2/3, 
and  set  /x  based  on  the  noise  level  in  the  observed  image  (/x  = 
104  for  all  results  reported  in  [1]). 


Fig.  1 .  Gradient  magnitude  “shrinkage”  function  for  a  =  2/3  and 
different  values  of  /3.  We  show  the  true  function,  as  well  as  the 
thresholded-linear  approximation  used  in  our  implementation. 


where  Wv(rx)  are  auxiliary  variables  constrained  to  lie  along 
VV(n).  The  solution  to  (8)  is  then  derived  using  an  iterative 
algorithm  that  alternately  updates  X  and  Wv(rx)  to  minimize 
C'(-),  starting  with  Xs  as  the  initial  estimate  for  X.  The 
parameter  is  increased  after  each  round  of  updates  to  Wv 
and  X,  by  a  factor  of  4  in  our  implementation.  We  ensure  that 
the  iterations  end  with  =  100 /x,  and  find  8  iterations  to  be 
sufficient  for  convergence. 


Updating  X:  Given  the  current  estimates  of  IFv(n),  X  is 
updated  to  minimize  the  first  two  terms  of  (9)  as 


i  fEv^[v]V[4i}N  +v/m{i}fny{i}] 

Evl^[V]|2+M//?Wi}]|2 


(10) 

Note  that  this  update  requires  only  one  forward  and  one 
inverse  Fourier  transform  per  channel  at  each  iteration,  since 
the  second  term  of  the  numerator  needs  to  be  only  computed 
once  for  an  input  image,  and  the  two  terms  of  the  denominator 
can  be  pre-computed  with  knowledge  of  the  kernels  and  the 
resolution  of  Y(n). 


Updating  Wv:  Correspondingly,  each  of  the  auxiliary  variables 
Wv(n)  can  be  independently  updated  based  on  the  current 
estimate  of  X(n),  to  minimize  the  latter  two  terms  of  (9)  and 
satisfy  the  color-direction  constraint  as: 


(Project)  r 

=  <Xv(n),Vv(n)>, 

(ii) 

(Shrink)  s 

—  argmin  \s\a  +  ^ \s  —  r|2, 
s  2 

(12) 

Wv(n) 

=  s  Vv(n). 

(13) 

Krishnan  and  Fergus  [4]  describe  an  analytical  solution  for  (12) 
when  a  =  2/3,  as  well  as  a  faster  look-up  table-based  approach. 
As  illustrated  in  Fig.  1,  we  find  that  the  relationship  between 
\s\  and  | r |  can  be  well  approximated  by  a  thresholded-linear 
function  (for  each  value  of  /?).  We  use  this  approximation  in 
our  implementation  since  it  is  faster  to  compute,  and  find  that 
it  gives  near  identical  results  as  the  analytical  solution. 


2.1  Half-quadratic  Splitting 

To  solve  the  optimization  problem  in  (8),  we  closely  follow  the 
approach  in  [4]  and  introduce  a  new  cost-function: 

C(X(n),  Wv(n),j9)=t  |  ^  [y{i}  (n)  -  (x{i}  *  k{i})(n)]2 

n,i 

+  f  ^||WMn)-Xv(n)||2 

V,n 

+  ^|<Wv(n),VV(n)>r,  (9) 

V,n 


3  Comparisons 

The  proposed  algorithm  is  used  to  generate  all  the  deblurred 
results  in  [1].  In  Fig.  2,  we  show  a  close-up  of  one  of  these 
results  as  well  as  a  comparison  to  the  deblurred  image  that 
would  have  been  recovered  using  the  baseline  greyscale  half¬ 
splitting  algorithm  [4],  with  the  same  value  of  a.  The  latter 
essentially  corresponds  to  deconvolving  each  color  channel 
independently,  and  while  the  two  approaches  have  nearly  iden¬ 
tical  running  times,  the  addition  of  color-gradient  constraints 
removes  spurious  chromatic  effects  caused  by  ringing  artifacts 
at  different  spatial  frequencies  in  different  color  channels. 
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Independent  per-channel  deconvolution  Deconvolution  with  color  constraints 


Fig.  2.  Comparison  of  deconvolution  with  and  without  color  constraints,  for  a  real  captured  image  from  [1].  We  see  that  adding 
color  constraints  removes  various  spurious  colors  (near  the  eyes  and  mouth)  that  appear  in  the  case  of  independent  per-channel 
deconvolution. 


PSNR:  23.76  dB  PSNR:  24.85  dB  PSNR:  24.93  dB  PSNR:  28.78  dB 


Fig.  3.  Comparison  to  results  reported  in  [5],  for  a  spectrally-uniform  kernel.  Despite  being  significantly  less  expensive 
computationally,  the  proposed  algorithm  yields  results  with  similar  quality  to  the  color  ILRS-based  method  of  Joshi  et  al.  [5]. 


We  next  show  results  on  images  from  the  Berkeley  seg¬ 
mentation  database  [6],  that  were  synthetically  blurred  with 
a  spectrally-uniform  kernel  and  used  for  evaluation  in  [5]. 
Figure  3  compares  the  deconvolution  results  from  our  method 
to  those  reported  for  the  color  IRLS-based  approach  of  [5], 
with  PSNR  values  for  both.  The  proposed  method  is  found  to 
have  equivalent  performance  to  [5],  while  offering  a  substantial 
advantage  in  terms  of  computational  efficiency 
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